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Abstract 



Extensive Monte-Carlo simulations were performed to study bond percola- 
tion on the simple cubic (s.c), face-centered cubic (f.c.c), and body-centered 
cubic (b.c.c.) lattices, using an epidemic kind of approach. These simulations 
provide very precise values of the critical thresholds for each of the lattices: 
p c (s.c.) = 0.248 812 6 ±0.000000 5, p c (f.c.c.) = 0.120 163 5 ± 0.000 001 0, and 
p c (b.c.c.) = 0.180 287 5 ± 0.0000010. For p close to p c , the results follow the 
expected finite-size and scaling behavior, with values for the Fisher exponent 
r (2.189 ± 0.002), the finite-size correction exponent Q (0.64 ± 0.02), and the 

scaling function exponent a (0.445 ± 0.01) confirmed to be universal. 
PACS numbers(s): 64.60Ak, 05.70. Jk 
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I. INTRODUCTION 



Percolation theory is used to describe a variety of natural physical processes, which have 
been discussed in detail by Stauffer and Aharony [[[]] and Sahimi [Q]. In two-dimensional 
percolation, either exact values or precise estimates are known for the critical thresholds 
and other related coefficients and exponents [H-flj. 

However, three-dimensional lattices are relevant for most natural processes. The most 
common of these are the simple cubic (s.c), the face-centered cubic (fx. a), and the body- 
centered cubic (b.c.c.) lattices. The percolation thresholds for these lattices are not known 
exactly, and the estimates that have been determined for the latter two lattices are much 
less precise than the values that have been found for typical two-dimensional systems. 

An exception is the case of the s.c. lattice, where fairly precise values have been deter- 
mined. A number of years ago, Ziff and Stell carried out a study of three-dimensional 
percolation for the site and bond percolation on the s.c. lattice [[7|. The values p c = 
0.248 812 ± 0.000 002 and p c = 0.311605 ± 0.000 010 were found for bond and site per- 
colation, respectively. The behavior for p away from p c was also studied, and it was found 
that the critical exponents can be described by the consistent set r = 116/53 ~ 2.1887, 
a = 24/53 « 0.4528, and u = 7/8, with errors ±0.001, ±0.001, and ±0.008, respectively. 
However, this work, which used a growth method essentially equivalent to the epidemic 
approach, remained unpublished as an internal report only |7j]. 

In a more recent work, Grassberger found p c = 0.248 814±0.000 003 and p c = 0.311 604± 
0.000 006 for bond and site percolation (respectively) on the s.c lattice ||, using an epidemic 
analysis |J. The precise agreement between these two independent works provides stong 
evidence that the procedures, random- number generators, error analyses, etc., implemented 
by both of these groups of authors are correct. Previous to these two works, p c was known 
only to about four significant figures, as discussed below. It is important to know p c very 
precisely so that other critical properties can be studied without bias. 

For bond percolation on the f.c.c. and b.c.c. lattices, the most accurate values appear to 
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be p c = 0.1200 ± 0.0002 and p c = 0.1802 ± 0.0002, respectively, recently found by van der 
Marck [[10] using the average crossing probability method |1[. Additional literature values 
are listed below. These thresholds are not precise enough for a study we have been carrying 
out on the universal excess cluster numbers in 3-D systems |LT] , and as a consequence we 
have conducted the present work. 

Here, we use a growth or epidemic analysis to obtain new, high-precision values for 
the percolation thresholds of the s.c, f.c.c, and b.c.c. lattices. We reaffirm that three- 
dimensional percolation follows the hypothesized finite-size and scaling behavior, and esti- 
mate the exponents and coefficients that enter into the finite-size and scaling functions. 

In the following sections, we discuss the simulation that we used to grow the percolation 
clusters and obtain our data. Then, we present and briefly discuss the results that we 
obtained from our simulations. 

II. SIMULATION METHOD 

We used a Monte-Carlo simulation of bond percolation on each of the three-dimensional 



lattices. This simulation employed the so-called Leath growth algorithm |L2[ to generate 
individual percolation clusters. The cluster was started at a seeded site that was centrally 
located on the lattice. From this site, a cluster was grown to neighboring sites by occupy- 
ing the connecting bonds with a certain probability p or leaving them unoccupied with a 
probability (1 — p). The unit vectors that were used to locate the neighboring sites for each 
lattice are summarized in Table I. Each of these clusters were allowed to grow until they 
reached an upper cut-off, at which point any cluster that was still growing was halted. This 
cut-off was 2 21 (2,097,152) wetted sites for the s.c. lattice and 2 20 (1,048,576) wetted sites 
for the f.c.c. and b.c.c. lattices. The size of the cluster was characterized by the number of 
wetted sites rather than the occupied bonds of the clusters, since the former figures directly 
in the growth algorithm procedure. 



In order to grow such large clusters, our simulation used the data-blocking scheme JIB], 
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in which large lattice is divided into smaller blocks and memory is assigned to a block 
only after the cluster has grown into it. In this case, the lattice, which has dimensions of 
2048 x 2048 x 2048, was divided into 64 3 blocks of dimensions of size 32 x 32 x 32. Bit 
mapping was also used to reduce the memory load of the large lattices. The upper six bits 
of each coordinate were used to tell where in the memory that block is mapped. The lower 
five bits of the y- and ^-coordinates tell the memory which word and the lower five bits of 
the x-coordinate tell the memory which bit of that word is used to store the site. With such 
a large lattice and the cut-offs we used, none of the clusters saw the lattice boundary, so 
there were no boundary effects. 

The simulation counted the number of clusters that closed in a range of (2 n , 2 n+1 — 1) sites 
for n = 0, 1, . . ., and recorded this number in the nth bin. If a cluster was still growing when 
it reached the upper cut-off, it was counted in the last (20th or 21st) bin. The simulation also 
kept track of the average number of occupied bonds and the average number of unoccupied 
bonds for clusters in each bin range, for reasons discussed below. 

Random numbers were generated by the four-tap shift-register rule x n = x n _47i © 
^n-i586 © a^n-6988 © ^n-9689 where © is the exclusive-or operation, which we have used (with 
apparent success) in numerous previous studies (e.g., [|§,[|). 



III. RESULTS 

A. Fisher exponent r 

Using the data obtained from the simulation, the number of clusters grown to a size 
greater than or equal to size s can be deduced. The probability of growing clusters with the 
number of wetted sites > s, P(s,p), when operating at the critical threshold (p = p c ) and 
neglecting finite-size effects, is predicted to follow [[[]: 

P(s, Pc )~As 2 -t (1) 
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In Fig. 1, data from our simulation plotted on a log- log plot shows good agreement with this 
equation. However, for small clusters, there is a clear nonlinear region, which represents the 
finite-size effect. For large s, the curves also become nonlinear, resulting from growing the 
clusters at values of p away from p c . Values of p that are greater than p c produce curves 
that increase for large s, and values of p below p c cause curves to decrease. 

Such behavior is similar to that seen in an epidemic analysis of a transition of interact- 
ing particle systems, and indeed, the cluster growth procedure is in fact essentially such a 
process. The only difference is that, in the usual epidemic analysis ||, one plots the growth 
probability as a function of the time (generation or chemical distance number for percola- 
tion), while here we plot that probability as a function of the size. But since the size scales 
with time, the two approaches are equivalent. 

The slope of the linear portion of the curves shown in Fig. 1 is equal to (2 — r), where 
t is the Fisher exponent JT4|. However, due to the nonlinear portions of the curves and the 
uncertainty of which p is precisely at the critical threshold, it is difficult to accurately deduce 
r from these plots. To make the behavior more pronounced, we plot difference between the 
data and a straight line of slope 2 — r in Fig. 2, for the f.c.c. lattice and different values of 
p. Here, the correct value of r produces a horizontal central portion of the curve. Applying 
this analysis for all three lattices yields r = 2.189 ± 0.002. Fig. 3 shows the effect of using 
slightly larger or smaller values of r for the f.c.c. lattice. 



In order to account for the small clusters, the finite-size correction must be added to (HI). 
Then P is predicted to follow 



B. Finite-size corrections. 




(2) 



where Q is the first correction-to-scaling exponent ||15|| . Like (|1|), (0) is only valid at the 
critical threshold. At the correct value of Q, (0) predicts there will be a linear relationship 
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between s T ~ 2 P(s,p c ) and s~ n . Fig. 4 shows plots between these two quantites for the s.c, 
f.c.c., and b.c.c. lattices, respectively. The values of Q which produced the best linear fit are 
summarized in Table II. The slope of the curves gives the value for the coefficient B in (0) 
and the intercept gives A in ([I]) and (fj). The values we found are summarized in Table III. 
Note that the last data point is shown on each of these plots, but they have been left out of 
linear fits because they represent clusters of only two sites, too small to be describe by just 
the first term of the finite-size corrections series. 

C. Percolation thresholds. 

In Fig. 1, we show the effect that values of p away from p c have on the data for large 
clusters. In order to account for the behavior when p ^ p c , a scaling function needs to be 
included in P. The behavior is then represented by 

P(s, Pc )~As 2 ~ T f((p-p c )s°) (3) 

which is valid rigorously in the scaling limit as s — > oo, p — > p c , such that (p — p c )s a = 
constant [IJ. Because p is close to p c , we can expand f(x) in a Taylor series: 

f((p ~ Pc)s a ) ~ 1 + C(p - p c )s° + .... (4) 

Eqs. (HD and (§]) describe how s T ~ 2 P(s,p) deviates from a constant value for large s when 
p is close to p c . Figs. 5(a), 5(b), and 5(c) show the plots of s T ~ 2 P(s,p) vs. s a for the s.c, 
f.c.c, and b.c.c. lattices, respectively. For these plots, we used the value of a = 0.453 from 
0, and the linearity of the plots confirms that it is a good value. These plots show a steep 
decline, which is the finite-size effect, followed by a mostly linear region as predicted by ([|). 
The linear portions of the curve become more nearly horizontal as the chosen value of p gets 
closer to p c . The value for the critical threshold can be deduced by interpolating from these 
plots. For the three lattices considered here, we thus find: 

p c (s.c.) = 0.248 812 6 ± 0.000 000 5 (5a) 
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p c (f.c.c.) = 0.120 1635 ± 0.000 001 (5b) 
p c (b.c.c.) = 0.180 287 5 ± 0.000 001 (5c) 

A more formal method of carrying out this interpolation can be obtained by plotting the 
slopes of the curves, shown in Figs. 5(a), 5(b), and 5(c), vs. p. From (^) and (^), it follows 
that: 

*^»= C <,-») + ... (6) 

which implies that the value of the critical threshold can be calculated from the x-intercept 
of this plot, as shown in Fig. 6. Figs. 5(a) and 6 yield consistent values of the critical 
threshold for the s.c. lattice. 



D. Scaling function. 

The linearized scaling function, as shown in (§) and (|J), can be studied efficiently by the 
following procedure. From (§) and it follows that: 

dP{ Q p P) = ACs 2 -^ + ... (7) 

Eq. (0) shows that we could determine 2 — r + a and AC directly from a measurement of 
the derivative dP/dp function of s. 

To develop a formula for this derivative, we consider the formal expression for the prob- 
ability of growing a cluster containing n occupied bonds and t vacant (perimeter) bonds: 

Pn^gn^il-PY (8) 

where g n ^ is the number of distinct clusters configurations with n occupied and t vacant 
bonds. From (ID, we can write the probability of growing a cluster of size greater than or 
equal to s as 

p(s, P ) = J2Y,9n,tP n (i-pY (9) 

n t 



where the sum is over all n and t such that the number of wetted sites is > s. 
Differentiating (0) with respect to p gives: 

8P(s,p) s t fn t 



n t \P 1 



dp nt' \P 

which can be simplified to the final form: 

dP(s,p) (n) (t) 



dp p 1 — p ' 

where (n) and (t) are the average number of occupied and unoccupied bonds, respectively, 
over all clusters of size > s. For a Monte-Carlo estimate of this derivative, we simply use (n) 
and (t) averaged over the sample of clusters. The values for these averages were recorded 



(10) 



(11) 



by our simulations, and were then used with (|TT|) to calculate the derivative. 



Fig. 7 shows a plot of the measured derivative (|TTD vs. s for the f.c.c. and b.c.c. lattices. 
From the slopes and intercepts, and using the values of A and r determined above, we can 
deduce a and C . For a we find 0.445 ± 0.01 for both cases, which is consistent with |7|], and 
the values of C we find are given in Table III. 



Another application of the derivative ( |TT| ) is that it can be used to produce new curves 
at values of p near p c that resemble the curves seen in Figs. 5(a), 5(b), and 5(c), by virtue 
of a Taylor expansion, 

dP(s,p) 



P(s, p 2 ) = P(s, pi) + (p 2 - Pi) 



dp 



(12) 



p=pi 

where the last derivative is determined using (|11|). This 'shifting' of the data is useful for 
correcting results taken close to, but not precisely at, p c , without carrying out a new set of 
time-consuming simulations. 



IV. DISCUSSION OF RESULTS 

A. Fisher exponent r 

Our results show excellent agreement with the expected relationship for the number of 
cluster of size greater than or equal to s, which is shown in (|l]). The only nonlinear portions 
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of the curve on a log-log plot occur when s is either small or large, which is where ([!]) is not 
valid. These nonlinearities are caused by the finite-size effect (small s) and the departure of 
p from p c (large s). 

The value of r was found by comparing curves that were similar to those shown in Figs. 
2 and 3, finding the value that produced the most nearly horizontal curve. As can be seen 
in Figs. 2 and 3 for the f.c.c. lattice, the reported value of r = 2.189 @ gives the best 
horizontal curve. This value of r also provides the best horizontal curves for the other two 
lattices, showing universality (with error bars of ±0.002). This value of r is also consistent 
with the values given by Stauffer and Aharony (2.18) [0J and Nakanishi and Stanley (2.19) 

B. Finite-size corrections. 

In Fig. 4, our data show excellent agreement to the prediction described by (|2[). It is 
expected that the exponent Q is a universal quantity while the coefficients A and B are 
different for the three different lattices. These expectations were born out by our results 
showing the best fit of data for all three lattices occurred at approximately the same value of 
Q (0.64 ±0.02), as shown in Table II, but different values of the coefficients as given in Table 
III. Our value for Q is significantly larger than the value reported by Nakanishi and Stanley 
(Q = 0.40) flBH - We believe this discrepancy is due to their using less precise values of the 
critical exponents, which leads to more inherent error in Q, and carrying out significantly 
fewer simulations, as were feasible at the time that that work was done. 

C. Percolation thresholds. 

In Figs. 5(a), 5(b), and 5(c), our data show excellent agreement with the predicted 
relationship as described by (|3[) and (^). The curves show the finite-size effects for small s 
and then become linear. The linear portion of the curves become more nearly horizontal as 
the value of p approaches p c , which is also predicted by the equations. 
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Our values for the critical thresholds ([5]) are consistent with most previous works. The 
critical threshold for the f.c.c. lattice has been reported as 0.1200 ± 0.0002 [[K]] and 0.119 



lyj. The critical threshold for the b.c.c. lattice has been reported as 0.1802 ± 0.0002 flC 
and 0.1803 [[[J. The critical threshold for the s.c. lattice has been reported as 0.248 812 ± 
0.000 002 0, 0.248 814 ± 0.000 003 §, 0.2488 §, 0.24875 ± 0.00013 [0, 0.2488 ± 0.0002 
and 0.2487 ± 0.0002 0. A number of years ago the higher values 0.2492 ± 0.0002 



and 0.2494 ± 0.0001 [^(J were reported, but these have since been recognized as probably 
erroneous due to a flaw in programming ||21|| . 



D. Scaling function. 

One expects that the value of the scaling function exponent a should be a universal 
quantity while the coefficient C should be different for each of the lattices. Values of these 
quantities were only calculated for the f.c.c. and b.c.c. lattices because in the s.c. lattice 
simulations we did not record the values of (n) and (t). The same value of a (0.445 ± 0.01) 
produced the best fit for both lattices. This value is significantly smaller than the value 
reported by Nakanishi and Stanley (0.504 ± 0.030) JOJ, but it is consistent with the values 
given by Ziff and Stell (0.453 ± 0.001) and Stauffer and Aharony (0.45). The values 
for the coefficient C which are summarized in Table III, were different for the two lattices. 
Thus, these results confirm the expectations. 



V. CONCLUSION 

Our work has produced the bond percolation critical threshold values given in (||). For 
the f.c.c. and b.c.c. lattices these results are at least two orders of magnitude more precise 
than previous values, while for the s.c. lattice the result is four times more precise. We 
find critical exponents consistent with those given in J!]], and a more precise value of the 
correction-to-scaling exponent Q than previously known. This increased precision is the 
result of having conducted extensive simulations with an inherently efficient procedure (the 
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epidemic method ||) along with programming techniques that allow one to simulate very 
large lattices ||13| . The results of this work allow other properties of three-dimensional 
percolation on these lattices to be studied equally precisely 

Because all of this work was performed in a relatively short amount of time, we believe 
that the epidemic approach is a more efficient way to find the critical threshold than the 
conventional crossing-probability methods. However, we did not make a direct time com- 
parison of two methods. Note that the epidemic growth algorithm used here can also be 
used to study systems of any dimension, unlike hull methods [p , 23fl , which are limited to 
two dimensions. 

Our measurements of P(s,p) for all s (except the smallest values) and for p very close 
to p c can be summarized in a single equation 

P(s,p) = s T ~ 2 (A + Bs- n )(l + C(p - p c )s°) (13) 

with all constants and exponents determined by our simulations. Although this equation 
mixes finite-size and bulk scaling forms, it provides a very accurate fit of the data in the 
regime we studied. 
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FIGURES 

FIG. 1. Plot of the raw simulation results on a log-log plot. This data is from the f.c.c. lattice 
and each of the five curves represent a different value of p: 0.1206 (square), 0.1204 (diamond), 
0.1202 (circle), 0.1200 (triangle), and 0.1198 (cross). Nonlinearities are caused by finite-size effects 
(small s) and p not equalling p c (large s). 

FIG. 2. Determination of the Fisher exponent, r, by plotting In s T ~ 2 P(s,p c ) vs. In s. These 
curves show data from the simulations of a f.c.c. lattice and each curve represents a different value 
ofp: 0.120165 (diamond), 0.120162 (circle), and 0.120 160 (triangle). These curves were produced 
using a value of r = 2.189, which produced the best horizontal curve for all three lattices. 

FIG. 3. Effect of varying r on the central curve of Figure 2 (f.c.c. lattice). These curves (all 
at p = 0.120 162, which is close to p c ) compare three different values of r: 2.192 (diamond), 2.189 
(circle), and 2.186 (square). The best horizontal fit is clearly produced by usingthe middle value. 

FIG. 4. Plot to determine the finite-size correction predicted by (2), for the s.c. (triangle), 
b.c.c. (circle), and f.c.c. (diamond) lattices. The values of $7, A, and B are summarized in Tables II 
and III. The last points of each curve, which represented clusters of just two wetted sites, were left 
out of the curve fit. Deviations from linearity for large r are caused by p being slightly different 
from p c . 

FIG. 5. Plot of the data vs. s CT . The linearity for large r demonstrates the validity of (3) 
and (4). The value of p that produces the best horizontal fit yields the critical threshold. Figure 
5(a) shows the data from the s.c. lattice, and each curve shows a different value of p: 0.248 820 
(diamond), 0.248 814 (square), 0.248 810 (circle), and 0.248 800 (triangle). Figure 5(b) shows the 
data from the f.c.c. lattice, and the values of p shown are: 0.120 165 (diamond), 0.120 162 (circle), 
and 0.120160 (triangle). Figure 5(c) shows the data from the b.c.c. lattices, and the values of p 
are: 0.180 300 (diamond), 0.180 2875 (circle), and 0.180 275 (triangle). 
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FIG. 6. Plot of the slope of the linear portion of the curves shown in Figure 5(a) as a function 
of p. The x-intercept gives the value of the critical threshold for the s.c. lattice. 

FIG. 7. Determination of a and C using dP/dp calculated from (11). By fitting the last fifteen 
data points to (7), 2 — r + a can be found from the slope of these curves and AC can be found 
from the y-intercept. The resulting values for a and C (using r and A determined above) for the 
f.c.c. (triangle) and b.c.c. (diamond) lattices are summarized in Tables II and III, respectively. 
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TABLES 



TABLE I. Unit vectors used to describe the neighbors in the s.c, f.c.c, and b.c.c. lattices, 
lattice vectors 

s.c. (1, 0, 0), (0, 1, 0), (0, 0, 1), (-1, 0, 0), (0, -1, 0), (0, 0, -1) 

f.c.c. (1, 1, 0), (1, -L 0), (-1, -L 0), (-1, 1, 0), (1, 0, 1), (-1,0, 1), 

(1, 0, -1), (-1,0, -1), (0, 1, 1), (0, -1, 1), (0, 1, -1), (0, -1, -1) 
b.c.c. (1, 1, 1), (1, 1, -1), (-1, 1, 1), (-1, 1, -1), (-1, -1, 1), (-1, -1, -1), (1, -1, 1), (1, -1, -1) 



TABLE II. Values of the universal finite-size correction exponent, and the scaling function 
exponent, a, for the s.c, f.c.c, and b.c.c. lattices, (a was not found for the s.c. lattice.) 
lattice Q a 

s.c. 0.63 

f.c.c. 0.66 0.4453 

b.c.c. 0.66 0.4433 



TABLE III. Non-universal coefficients for finite-size correction and the scaling function of the 
s.c, f.c.c, and b.c.c. lattices. (C was not found for s.c. lattice.) 

lattice ABC 
s.c. 0.813 0.178 

f.c.c. 0.767 0.186 7.95 

b.c.c. 0.776 0.187 5.57 
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Figure 1 - Lorenz, Ziff 




Figure 2 - Lorenz, Ziff 




Figure 3 - Lorenz, Ziff 
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Figure 4 - Lorenz, Ziff 




Figure 5a - Lorenz, Ziff 




Figure 5b - Lorenz, Ziff 




Figure 5c - Lorenz, Ziff 




Figure 6 - Lorenz, Ziff 
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Figure 7 - Lorenz, Ziff 




